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The exact stochastic decomposition of non-Markovian dissipative quantum dynamics is combined 
with the time-dependent semiclassical initial value formalism. It is shown that even in the challeng- 
ing regime of moderate friction and low temperatures, where non-Markovian effects are substantial, 
this approach allows for the accurate description of dissipative dynamics in anharmonic potentials 
over many oscillation periods until thermalization is reached. The problem of convergence of the 
stochastic average at long times, which plagues full quantum mechanical implementations, is avoided 
through a joint sampling of the stochastic noise and the semiclassical phase space distribution. 
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The understanding of the nonequilibrium dynamics of 
open quantum systems has been a central challenge in the 
last decades [HISI- In recent years the subject has gained 
considerable interest due to experimental progress which 
allows for the tailoring and manipulation of quantum 
matter on ever larger scales. In mesoscopic physics, for 
instance, superconducting circuits have been realized to 
observe coherent dynamics and entanglement Q . Similar 
advance has been achieved on molecular scales with the 
detection of interferences in wave packet dynamics and 
the control of the population of specific molecular states 
Typically, these systems are in contact with a large 
number of environmental degrees of freedom, e. g. elec- 
tromagnetic modes of the circuitry or residual vibronic 
modes, which give rise to decoherence and relaxation Q. 
In many cases, the idealization of an isolated system must 
inevitably be replaced by an open-system theory. 

In the standard approach to open quantum systems, 
the reduced dynamics of the system of interest is ob- 
tained by tracing out "reservoir" degrees of freedom from 
the conservative system-plus-reservoir dynamics, e. g., 
through projection operator techniques. Alternatively, 
this program can be carried out through exact path inte- 
gral expressions for the reduced density matrix [6| , which 
became widely used in the 1980s [l|. The distinguish- 
ing feature of dissipative path integrals is an influence 
functional which describes self-interactions non-local in 
time. Hence, a simple quantum mechanical analogue to 
the classical Langevin equation is not known; commonly 
used equations, such as Master/Redfield equations [2| 
in the weak-coupling case and quantum Smoluchowski 
equations 0] for reservoir-dominated dynamics, rely on 
perturbation theory. In intermediate domains, quantum 
Monte Carlo techniques have been put forward for tight- 
binding systems, but the achievable propagation times 
are severely limited by the dynamical sign problem. 

Recently, it has been shown that the influence func- 
tional can be exactly reproduced through stochastic av- 
eraging of a process without explicit memory [1,0]. The 
formulation turned out to be particularly efficient for 
weak to moderate friction and low temperatures [3, [lo| , 



a regime which lies beyond the strict validity of Redfield 
equations on the one hand and beyond the applicability 
of Monte Carlo schemes on the other hand. The draw- 
back is though that for nonlinear systems the convergence 
of the stochastic average for relatively long times is still 
an unsolved problem. Some progress has been made for 
the spin-boson model by u sing a hierarchical approach 
to quantum memory terms [ll|. A reliable and efficient, 
generally applicable method to tackle the dissipative dy- 
namics for continuous systems in this challenging param- 
eter regime is still missing. 

Here we address this issue by combining the exact 
stochastic Schrodinger formulation with a semiclassical 
real-time technique based on a time-dependent initial 
value representation of the quantum mechanical prop- 
agator. The central finding is that this procedure cir- 
cumvents the main obstacle of the exact formulation and 
allows for accurate simulations up to times where equi- 
libration sets in. It is important to note that a direct 
stationary-phase evaluation of the double path integral 
for the reduced density is not a consistent semiclassical 
approximation since the classical limit of open-system dy- 
namics, the Langevin equation, is not recovered [ll IT^. 

Undamped systems allow a powerful method based 
on the semiclassical propagator of Herman and Kluk 
(HK) [isjl , which has seen an impressive number of appli- 
cations ranging from atomic [w| to chemical physics 15 1 
after the work of Kay [l^ stimulated renewed interest in 
the approach. The HK propagator was recently shown 
to be the leading term of consistent series expansions of 
the exact quantum propagator of a Hamiltonian system 



171 [18| . Currently, efforts are under way to directly ex- 
tend this approach to dynamics with quantum memory 
effects [31 . In this Letter, we will use a memory- free rep- 
resentation, which accounts for non-Markovian dynamics 
entirely through correlations of complex noise forces and 
thus suits numerical applications. 

We start with the standard decomposition of the total 
Hamiltonian 
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as a sum of a system part, that for reasons of simplicity 
shall here depend on one degree of freedom x, a bath part 
consisting of an infinity of harmonic oscillators together 
with a bilinear interaction between them. In case of a 
factorized initial density with a bath residing in thermal 
equilibrium at temperature T one derives a path integral 
expression for the time-evolved reduced density matrix 
of the form 

p{xf,x'j:,t) = dxidXip{xi,x'^,0) 'D[xi]'D[x2] 



exp |-(5s[xi] - Ss[x2])j F[xuX2U2) 

where the two real time paths xi and X2 run in time 
t from Xi and x'^ to Xf and x^, respectively. They are 
coupled by the influence functional, which takes the form 
F[y,r] — exp(— <!>[?/, r]/7i) with 

1 /■* r 

$[y,r] — — du dvy{u)[L' {u — v)y{v) + 2iL" {u — v)r 
ri- Jo Jo 

/ duy{u)r{u) , 

where y = xi — X2^ r = (xi + X2)/2 denote differ- 
ence and sum paths, respectively. The complex valued 
friction kernel L{t) = L'{t) -\- iL"{t) is related to the 
force-force auto-correlation function of the bath and com- 
pletely determined by its spectral density J{uj) and in- 
verse temperature /3. The static susceptibility denoted by 
/I = — duL" (u) / {2h) is a property of the reservoir. 

In [9] it has been shown that a stochastic unraveling 
of the forward and the backward paths leads to 

p{xf,x'f,t) = Jdxi J dx'ip{xi,Xi,0) 

xM[K,,ixf,t; X,, G){K,,{x'f,t- x^Q))*] , (A) 

where M denotes the average over noise realizations Zj 
(j = l,2), with the noise augmenting the system actions 
via 



where f (i) = \[zi{t) + z*{t)] and v{t) = ^[zi{t) - z^{t)]. 
The reduced density matrix ([2]) is reproduced exactly 
by averaging p obtained from equations ([6]) and ([7|) 
when the correlations of ^ and v reproduce the inte- 
gral kernel of the influence functional: M[^{t)^{t')] — 
L'{t - t'), M[^{t)v{t')] = {2i/h)Q{t - t')L"{t - t'), and 
M[v{t)v{t')] = (8 denotes the Heaviside step function). 

The linear equations ([6]) and ([7|) are formally exact in 
the sense that their Monte Carlo sampling will eventually 
converge to yield Eq. ([2]). However, they are of limited 
use for practical calculations since individual samples do 
not stay normalized, which slows down convergence and 
makes a direct samping impractical This slowdown 

can be reduced [lO| by an exact mapping of the stochas- 
tic processes dU and ([7]) to a trace-conserving process 
(similar to a Girsanov transform [2l!|). Observing that 
only the last terms in the square brackets of equations 
([6]) and ([7]) lead to a change of tr /5, the first step of the 
transform consists of subtracting an arbitrary "reference 
('^)]trajectory" f„ from x in these terms jl^l- The effect of 
this modification is canceled exactly by a corresponding 
(|3^change in the probability measure, which can be repre- 
sented by the substitution 



i = i- I dux{t - u)ru 



Szj[xj\ = Ss[x.j\+ii I duxj{u)^+ duxj{u)zj{u) 
Jo Jo 



(5) 



in the path integral expressions of the respective propaga- 
tors ■ This stochastic unraveling differs from a similar 
one by Strunz et al. 20] through the appearance of two 
noise variables, allowing for the elimination of quantum 
memory effects. 

Representing a general initial density operator through 
p{t = 0) = |4'i)(5'2| (or through an ensemble of such 
projectors) leads to two Schrodinger equations 



(8) 



where xC"*^) = —Q{u)L"{u)/2T% is the response function 
of the reservoir. Details of the transform are given in 
Refs. 3 and 0. With f„ = (*i|x|^'2)„, the diffusion 
of tr p is eliminated. However, this can lead to subtle 
mathematical difficulties limiting the times for which nu- 
merical simulations are stable [10]. There are two situa- 
tions known to be free of such instabilities, namely, linear 
systems and the classical limit. The idea is thus to com- 
bine the stochastic quantum dynamics with a semiclas- 
sical propagation scheme based on the frozen Gaussian 
approximation pioneered by Herman and Kluk [l3| and 
Heher fisj. 

The propagation of individual samples of the stochastic 
processes ([6]) and ([7]) by the HK propagator differs from 
the evolution of a closed system only through the addi- 
tion of simple potential terms, up to quadratic order, to 
the system Hamiltonian. Hence the asymptotic conver- 
gence properties of the HK propagator for closed systems 
lis., 2A, 2S\ are 'inherited' by our stochastic samples. 

The semiclassical HK propagator is given in terms of 
a phase space integral as 



K{xf,t,Xi,0) 



dpidqi 
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{xf\9',{Pt,qt)) 



in\i>i) = 



Hs-C{t)x + ^x^ + K*{t)x 



*i> (6) 
*2), (7) 



X i?(p„ g„ t)e'^fe •''-*)/^ (g^ (p„ ) I a;.) . (9) 

Here complex valued Gaussian wave packets {x\g~f) ~ 
exp{— ^(a; — qY + yip{x — q)} of fixed width parame- 
ter 7 have been introduced, centered around the initial 
phase space points pi, qi and the time-evolved phase space 
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points Pt,qt, respectively. The pre-exponential factor R 
contains a complex valued combination of stability ma- 
trix elements and the action reads as in ([5]) with the 
replacement of the noise force described in equation ([S]). 
Obtaining the final density matrix p{xf,x'pt) involves 
three Monte Carlo integrations, two over the forward and 
backward phase spaces of the semiclassical propagators 
and an additional one over the noise trajectory distribu- 
tion. 

The classical trajectory entering equation ^ is ob- 
tained from the quasiclassical dynamics of the fixed- 
width Gaussians under the transformed versions of equa- 
tions (O and (O. The complex forces ^ and u do not 
extend the phase space to complex numbers: The frozen 
Gaussians are, up to a trivial phase factor, coherent 
states la) = e^l^'l/^ e"°* |0) with 



(10) 



It is clear from this equation that complex values of q and 
p lead only to states already described by a real-valued 
phase space. In compact form, the classical equations of 
motion derived from ^ and ([7]) read (j=l,2) 
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(11) 



with fi = ^ + and /2 = I* — f^^*, to be solved for 
real qj and Pj. Taking the limit h ^ 0, and integrating 
by parts in equation ([5]) , the classical Langevin equation 
is indeed recovered from equation (jlip . 

In the semiclassical context, the reference trajectory 
is again obtained by demanding that the z^-dependent 
terms in equations ([6]) and ([7]) do not change the trace 
of the sample, which leads to the simple condition 



(12) 



(ai+a^)/V27 



This definition of f„ makes reference to a single pair of 
semiclassical trajectories. It is therefore advantageous to 
merge the integrations over the two HK phase spaces and 
the function space of noise trajectories ^(i) and ^{t) into 
a joint Monte Carlo sampling scheme where new initial 
phase space points and a new noise trajectory are drawn 
for each sample (26j . 

To demonstrate the approach, in the remainder we will 
concentrate on the dynamics in a one-dimensional Morse 
potential 



V{x) = hD {1 — cxp[— a(a; — xo)/xo]}'^ 



(13) 



which may serve e. g. as a simple model for the rela- 
tive motion of a diatomic molecule with equilibrium sep- 
aration xo- The HK propagation for the isolated sys- 
tem is known to accurately reproduce the exact quantum 
dynamics. Here we introduce the scaled displacement 
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FIG. 1: Expectation value of position and variance in a dissi- 
pative Morse oscillator for decreasing temperature h/S = 0.5 
(full), hp — 1 (dotted), hl3 = 10 (dashed) at constant rj = 0.1. 
Arrows indicate thermal expectation values of the unper- 
turbed oscillator (see text). 



X = (x — xq)/xq and the frequency unit h/ {uixq). In the 
following we use the parameters D = 30 and a = 0.08 
corresponding to a total number of 97 bound states. In 
these units the frequency loq for small oscillations around 
the well minimum located at a; = is cjq = a\JlD = 0.62. 
In the Markovian case and with the restriction to the 
bound state part of the spectrum a similar system has 
been studied in [23]. As the initial state of the propa- 
gation we have chosen a Gaussian wave packet shifted 
away from the potential minimum with (i)o = 1 and 
zero initial momentum. The spectral density of the bath 
oscillators is taken to be Ohmic with an algebraic cutoff 



J(c.) = 



(l+^V^c'F 



(14) 



where the cutoff frequency lOc = 10 is well beyond any rel- 
evant system frequencies and the dimensionless coupling 
strength is denoted by 77. 

Let us first look at the temperature dependence of 
mean and variance in position for a fixed interaction 
strength r; = 0.1, displayed in Fig. [TJ As expected, 
damped oscillations occur for intermediate times. In 
sharp contrast to the full quantum mechanical case, how- 
ever, the simulations are stable over long periods of time 
until equilibration is approached. Arrows indicate equi- 
librium values obtained by Boltzmann weighting analytic 
bound-state results from [2^. At higher temperatures 
the mean saturates at a value considerably away from 
zero thus revealing the substantial influence of the an- 
harmonicity. This influence gradually decreases when 
the temperature is lowered. Only at the lowest tempera- 
ture (?i/3 = 10, dashed arrows) are the results consistent 
with a harmonic approximation of the oscillator. Sim- 
ilar effects are seen in the position variance, with the 
additional feature that coherent oscillations are strongly 
smeared out at higher temperatures. Importantly, for the 
lowest temperatures considered here, one has tOchP 3> 1 
so that retardation effects are strong and the dynamics 
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FIG. 2: Same as in Fig. [T] but for increasing coupling tj — 
(dashed), r] = 0.05 (dotted), 77 = 0.1 (full), all at hp = 1. 

is far from being Markovian. Further, since r]hp w 1 
we are outside the weak coupling regime in the sense of 
Master/Redfield equations. To check the thermal nature 
of the final expectation values, we used different initial 
preparations and found that the results converge to the 
same long time average (not shown) . The total number of 
sampling points needed to generate all presented results 
is 10^; the calculations take a few hours on a desktop PC. 

In Figure [2] the same expectation values are depicted 
at a fixed inverse temperature h/3 = 1 but varying dissi- 
pation strengths. Of course, pronounced coherent oscil- 
lations are seen for vanishing coupling, where again the 
anharmonicity of the potential shows up in smaller oscil- 
lation amplitudes for negative than for positive values of 
position. Now a marked difference between the isolated 
system and the case of weak damping can be observed: 
Even for weak damping the amplitudes of the oscilla- 
tions decrease strongly, with centers shifted towards the 
soft side of the potential. The variance increases as a 
function of time, while its initial oscillations decay on 
a somewhat shorter timescale. We mention that at the 
fixed temperature 7i/3 = 1, almost identical long time av- 
erages are reached for all small friction strengths that we 
have investigated. 

To summarize, we have introduced a well-defined semi- 
classical propagator for non-Markovian open quantum 
systems. It was derived by consistently applying a known 
stochastic decomposition of the influence functional ap- 
pearing in exact open-system path integrals. A test us- 
ing a standard model of molecular physics shows that the 
dissipative semiclassical propagator accurately describes 
the time evolution over long periods of time, up to ther- 
mal equilibration between open system and environment. 
With modest computational resources, we provide data 
for this model which are difficult if not impossible to 
obtain by other means in the low-temperature regime of 
non-Markovian dynamics. For systems where corrections 
to the HK propagator are needed, extensions according 
to [3| are straightforward. We estimate that the num- 
ber of samples needed for a higher-dimensional quantum 



system will not increase substantially beyond the num- 
ber of samples in the one-dimensional case because the 
phase-space sampling converges more rapidly than the 
noise average. 
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